Packages
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr")
if (!require("lme4")) install.packages("lme4")
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # p-values
if (!require("car")) install.packages("car")
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown formatBackground
This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis) (published in the Journal __ in 2023). Data was collected and analyzed by Calvin Davis and Savannah Weaver, under the advising of Dr. Emily Taylor at California Polytechnic State University.
Load Data
Get the RDS files created in the data wrangling Rmd:
temp_time_series <- read_rds("Data/temp_time_series.RDS")
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS")Stats & Models
Temp Correlation
Check the correlation between dorsal and cloacal temps:
temp_corr <- lm(data = temp_time_series,
calibrated_dorsal_temp ~ calibrated_cloacal_temp)
summary(temp_corr)##
## Call:
## lm(formula = calibrated_dorsal_temp ~ calibrated_cloacal_temp,
## data = temp_time_series)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0389 -0.7784 0.0075 0.5911 5.3351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.966310 0.049820 139.8 <2e-16 ***
## calibrated_cloacal_temp 0.747869 0.001899 393.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.388 on 25249 degrees of freedom
## Multiple R-squared: 0.86, Adjusted R-squared: 0.8599
## F-statistic: 1.55e+05 on 1 and 25249 DF, p-value: < 2.2e-16They are actually pretty different, so run parallel models of CEWL, one with dorsal temp, and one with cloacal temp.
CEWL (raw) ~ temp LMM
First, test linear effects of temperature:
CEWL_LMM_02a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
calibrated_cloacal_temp *
treatment +
(1|date/lizard_ID))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesCEWL_LMM_02b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
calibrated_dorsal_temp *
treatment +
(1|date/lizard_ID))Compare to the inclusion of polynomial factors:
CEWL_LMM_03a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_cloacal_temp, 2)*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_03a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02a 9 87929 88002 -43956 87911
## CEWL_LMM_03a 12 81676 81774 -40826 81652 6259.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_03b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_dorsal_temp, 2)*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_03b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02b 9 88775 88848 -44379 88757
## CEWL_LMM_03b 12 83092 83190 -41534 83068 5689.1 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The poly models have significantly better explanatory power than the linear models.
Also try log transformation to make sure a polynomial effect is best:
CEWL_LMM_04a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
(log(calibrated_cloacal_temp))*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_04a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 88930 89003 -44456 88912
## CEWL_LMM_02a 9 87929 88002 -43956 87911 1000.4 0anova(CEWL_LMM_04a, CEWL_LMM_03a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 88930 89003 -44456 88912
## CEWL_LMM_03a 12 81676 81774 -40826 81652 7259.7 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_04b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
log(calibrated_dorsal_temp)*treatment +
treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_04b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 89762 89835 -44872 89744
## CEWL_LMM_02b 9 88775 88848 -44379 88757 987.18 0anova(CEWL_LMM_04b, CEWL_LMM_03b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 89762 89835 -44872 89744
## CEWL_LMM_03b 12 83092 83190 -41534 83068 6676.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The logarithmic models are NOT better than the linear or polynomial ones.
Check assumptions:
plot(CEWL_LMM_03a)
plot(CEWL_LMM_03b)
vif(CEWL_LMM_03a)## GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_cloacal_temp, 2) 9.312576e+05 2 31.064721
## treatment 1.100911e+00 2 1.024326
## poly(calibrated_cloacal_temp, 2):treatment 9.312225e+05 4 5.573547vif(CEWL_LMM_03b)## GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_dorsal_temp, 2) 5.528187e+05 2 27.267523
## treatment 1.107683e+00 2 1.025897
## poly(calibrated_dorsal_temp, 2):treatment 5.527949e+05 4 5.221803They both have a really weird spike, but otherwise there’s no fanning or different pattern. VIF is negligible.
Re-run with p-values:
CEWL_LMM_besta <- lmerTest::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_cloacal_temp, 2) * treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_besta)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 81611.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6892 -0.6401 -0.0881 0.5378 7.8576
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 11.773 3.431
## date (Intercept) 23.794 4.878
## Residual 1.468 1.212
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 0.6557 2.4503
## poly(calibrated_cloacal_temp, 2)1 1422.2807 62.9113
## poly(calibrated_cloacal_temp, 2)2 -3419.1801 82.4963
## treatmentCooling 6.6956 1.5010
## treatmentHeating 16.7664 1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238 62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 3565.6587 82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577 62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 3338.2055 82.5256
## df t value Pr(>|t|)
## (Intercept) 5.4550 0.268 0.798840
## poly(calibrated_cloacal_temp, 2)1 25235.9702 22.608 < 2e-16
## poly(calibrated_cloacal_temp, 2)2 25238.1540 -41.446 < 2e-16
## treatmentCooling 30.3218 4.461 0.000104
## treatmentHeating 29.8965 10.547 1.36e-11
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling 25235.9622 -18.460 < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 25238.1521 43.205 < 2e-16
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating 25235.9530 -19.438 < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 25238.1526 40.451 < 2e-16
##
## (Intercept)
## poly(calibrated_cloacal_temp, 2)1 ***
## poly(calibrated_cloacal_temp, 2)2 ***
## treatmentCooling ***
## treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159
## ply(c__,2)2 0.164 -0.914
## tretmntClng -0.336 0.260 -0.268
## tretmntHtng -0.324 0.244 -0.251 0.505
## pl(__,2)1:C 0.159 -0.999 0.914 -0.259 -0.244
## pl(__,2)2:C -0.164 0.914 -1.000 0.268 0.251 -0.913
## pl(__,2)1:H 0.159 -0.999 0.914 -0.259 -0.244 0.999 -0.913
## pl(__,2)2:H -0.164 0.914 -1.000 0.268 0.251 -0.914 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.914#anova(CEWL_LMM_besta, type = "3", ddf = "Kenward-Roger")
CEWL_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_dorsal_temp, 2) * treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_bestb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 83028.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7628 -0.6345 -0.0833 0.5160 8.3769
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.996 3.162
## date (Intercept) 21.953 4.685
## Residual 1.554 1.246
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 3.248 2.336
## poly(calibrated_dorsal_temp, 2)1 1534.609 67.430
## poly(calibrated_dorsal_temp, 2)2 -2711.185 75.106
## treatmentCooling 3.529 1.386
## treatmentHeating 14.341 1.468
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765 67.457
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 2852.577 75.134
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961 67.469
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 2666.790 75.141
## df t value Pr(>|t|)
## (Intercept) 5.353 1.390 0.2194
## poly(calibrated_dorsal_temp, 2)1 25237.527 22.759 < 2e-16
## poly(calibrated_dorsal_temp, 2)2 25236.097 -36.098 < 2e-16
## treatmentCooling 30.609 2.546 0.0162
## treatmentHeating 30.174 9.766 7.46e-11
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling 25237.523 -20.024 < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 25236.096 37.967 < 2e-16
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating 25237.517 -20.246 < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 25236.103 35.490 < 2e-16
##
## (Intercept)
## poly(calibrated_dorsal_temp, 2)1 ***
## poly(calibrated_dorsal_temp, 2)2 ***
## treatmentCooling *
## treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162
## ply(c__,2)2 0.163 -0.946
## tretmntClng -0.327 0.273 -0.276
## tretmntHtng -0.315 0.257 -0.260 0.507
## pl(__,2)1:C 0.162 -1.000 0.946 -0.273 -0.257
## pl(__,2)2:C -0.163 0.946 -1.000 0.276 0.259 -0.946
## pl(__,2)1:H 0.162 -0.999 0.946 -0.273 -0.257 0.999 -0.945
## pl(__,2)2:H -0.163 0.946 -1.000 0.276 0.259 -0.945 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.945#anova(CEWL_LMM_bestb, type = "3", ddf = "Kenward-Roger")
# double check sample sizes
temp_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 11
## 2 Cooling 12
## 3 Heating 10These are the best models for CEWL ~ temperature.
Export:
broom.mixed::tidy(CEWL_LMM_besta) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")
broom.mixed::tidy(CEWL_LMM_bestb) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")CEWL (raw) ~ time LMM
NOTE: Start times vary from 97-170 seconds after starting time series.
Make a polynomial model first:
rates_LMM_01 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
drop1(rates_LMM_01)## Single term deletions
##
## Model:
## CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC
## <none> 96906
## poly(time_min, 2):treatment 4 111987Also make a linear model version:
rates_LMM_02 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
time_min * treatment +
(1|date/lizard_ID))
anova(rates_LMM_02, rates_LMM_01)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_02 9 100913 100987 -50447 100895
## rates_LMM_01 12 96906 97005 -48441 96882 4012.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Polynomial is much better than linear.
Also compare to log version:
rates_LMM_03 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
log(time_min) * treatment +
(1|date/lizard_ID))
anova(rates_LMM_02, rates_LMM_03)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_02 9 100913 100987 -50447 100895
## rates_LMM_03 9 97121 97196 -48552 97103 3791.4 0anova(rates_LMM_01, rates_LMM_03)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_03 9 97121 97196 -48552 97103
## rates_LMM_01 12 96906 97005 -48441 96882 220.96 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The polynomial model is best (better than both log and lm options).
Check linear regression assumptions:
plot(rates_LMM_01)
vif(rates_LMM_01)## GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2) 9.050897 2 1.734494
## treatment 1.000002 2 1.000000
## poly(time_min, 2):treatment 9.050904 4 1.317002Again, there is a really weird spike, but otherwise there’s no fanning or pattern.
re-run with p-values:
rates_LMM_best <- lmerTest::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_best)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## Data: CEWL_time_series
##
## REML criterion at convergence: 96853.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3325 -0.5700 -0.0438 0.5178 8.0025
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.976 3.158
## date (Intercept) 17.817 4.221
## Residual 1.754 1.324
## Number of obs: 28405, groups: lizard_ID:date, 37; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 14.957 2.099 5.189 7.126
## poly(time_min, 2)1 -137.974 2.309 28362.024 -59.765
## poly(time_min, 2)2 51.100 2.322 28362.016 22.008
## treatmentCooling -8.473 1.297 30.005 -6.533
## treatmentHeating 3.416 1.280 30.061 2.669
## poly(time_min, 2)1:treatmentCooling -88.851 3.376 28362.571 -26.319
## poly(time_min, 2)2:treatmentCooling 83.367 3.337 28362.083 24.980
## poly(time_min, 2)1:treatmentHeating 301.172 3.201 28362.093 94.092
## poly(time_min, 2)2:treatmentHeating -108.483 3.201 28362.048 -33.893
## Pr(>|t|)
## (Intercept) 0.000723 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 < 2e-16 ***
## treatmentCooling 3.17e-07 ***
## treatmentHeating 0.012151 *
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 0.000
## ply(tm_,2)2 0.000 -0.067
## tretmntClng -0.305 0.001 0.000
## tretmntHtng -0.317 0.001 0.000 0.499
## ply(_,2)1:C 0.000 -0.684 0.046 0.000 0.000
## ply(_,2)2:C 0.000 0.047 -0.696 0.000 0.000 0.008
## ply(_,2)1:H 0.000 -0.721 0.049 0.000 0.000 0.493 -0.034
## ply(_,2)2:H 0.000 0.049 -0.725 0.000 0.000 -0.033 0.505
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.029#anova(rates_LMM_best, type = "3", ddf = "Kenward-Roger")
# double check sample sizes
CEWL_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 12
## 2 Cooling 12
## 3 Heating 13Export:
broom.mixed::tidy(rates_LMM_best) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_time.csv")Body Temperature ~ Time
linear models:
body_temp_LMM_01a <- lme4::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
time_min * treatment +
(1|date/lizard_ID))
body_temp_LMM_01b <- lme4::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
time_min * treatment +
(1|date/lizard_ID))polynomial models:
body_temp_LMM_02a <- lme4::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
anova(body_temp_LMM_02a, body_temp_LMM_01a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## body_temp_LMM_01a: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02a: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## body_temp_LMM_01a 9 63638 63712 -31810 63620
## body_temp_LMM_02a 12 41291 41388 -20633 41267 22354 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1body_temp_LMM_02b <- lme4::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
anova(body_temp_LMM_02b, body_temp_LMM_01b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## body_temp_LMM_01b: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02b: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## body_temp_LMM_01b 9 77531 77604 -38756 77513
## body_temp_LMM_02b 12 71434 71532 -35705 71410 6102.6 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Again, the polynomial models are MUCH better.
body_temp_LMM_besta <- lmerTest::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_besta)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 41254.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8989 -0.4649 -0.0011 0.3942 6.3225
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 3.0570 1.7484
## date (Intercept) 0.1856 0.4308
## Residual 0.2967 0.5447
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 27.4181 0.5621 16.5563 48.780
## poly(time_min, 2)1 52.0361 0.9376 25212.0138 55.500
## poly(time_min, 2)2 6.0643 0.9435 25212.0089 6.428
## treatmentCooling -4.0686 0.7317 26.4065 -5.560
## treatmentHeating -0.7262 0.7697 27.6621 -0.944
## poly(time_min, 2)1:treatmentCooling -764.1781 1.3382 25212.2800 -571.054
## poly(time_min, 2)2:treatmentCooling 161.4913 1.3255 25212.0407 121.838
## poly(time_min, 2)1:treatmentHeating 716.1235 1.3557 25212.0589 528.233
## poly(time_min, 2)2:treatmentHeating -63.6826 1.3531 25212.0326 -47.063
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 1.32e-10 ***
## treatmentCooling 7.32e-06 ***
## treatmentHeating 0.354
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.676 0.000 0.000
## tretmntHtng -0.647 0.000 0.000 0.490
## ply(_,2)1:C 0.000 -0.701 0.054 0.000 0.000
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.000 -0.692 0.053 0.000 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021write.csv(broom.mixed::tidy(body_temp_LMM_besta),
"./Results_Statistics/temp_time_best_model_clotemp.csv")
body_temp_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_bestb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 71392.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8879 -0.3478 -0.0225 0.3300 6.0970
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 1.5384 1.2403
## date (Intercept) 0.3563 0.5969
## Residual 0.9812 0.9906
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 27.7236 0.4606 11.4180 60.191
## poly(time_min, 2)1 48.0115 1.7051 25212.0897 28.158
## poly(time_min, 2)2 0.6860 1.7158 25212.0585 0.400
## treatmentCooling -2.8157 0.5208 26.3868 -5.406
## treatmentHeating -1.5969 0.5509 27.0887 -2.899
## poly(time_min, 2)1:treatmentCooling -697.0769 2.4336 25213.7698 -286.444
## poly(time_min, 2)2:treatmentCooling 133.9832 2.4104 25212.2645 55.585
## poly(time_min, 2)1:treatmentHeating 484.8482 2.4654 25212.3649 196.660
## poly(time_min, 2)2:treatmentHeating -42.9964 2.4608 25212.2024 -17.473
## Pr(>|t|)
## (Intercept) 1.19e-15 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 0.68930
## treatmentCooling 1.10e-05 ***
## treatmentHeating 0.00734 **
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.584 0.001 0.000
## tretmntHtng -0.560 0.001 0.000 0.480
## ply(_,2)1:C 0.001 -0.701 0.054 0.000 -0.001
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.001 -0.692 0.053 -0.001 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021write.csv(broom.mixed::tidy(body_temp_LMM_bestb),
"./Results_Statistics/temp_time_best_model_dorstemp.csv")Rate of Change
We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.
First, calculate rate of change for each lizard:
change_LM <- lm(data = CEWL_time_series,
CEWL_g_m2h ~
time_min*lizard_ID)Now get the trend/slope for each individual and make into a variable in df.
lizard_coeffs <- data.frame(emtrends(change_LM, "lizard_ID", var = "time_min")) %>%
left_join(read_rds("Data/collection_dat_formatted.RDS"),
by = c("lizard_ID" = "individual_ID")) %>%
dplyr::select(lizard_ID, rate_change = time_min.trend,
SE, df, lower.CL, upper.CL,
date,
treatment,
mass_g, SVL_mm,
chamber_time_elapsed_hr
) %>%
dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs)## lizard_ID rate_change SE df
## 401 : 1 Min. :-0.8023 Min. :0.008416 Min. :28331
## 402 : 1 1st Qu.:-0.3728 1st Qu.:0.008672 1st Qu.:28331
## 404 : 1 Median :-0.1621 Median :0.008739 Median :28331
## 406 : 1 Mean :-0.1074 Mean :0.009524 Mean :28331
## 407 : 1 3rd Qu.: 0.1088 3rd Qu.:0.008920 3rd Qu.:28331
## 408 : 1 Max. : 0.6960 Max. :0.024080 Max. :28331
## (Other):31
## lower.CL upper.CL date treatment
## Min. :-0.81947 Min. :-0.78520 Min. :2021-10-05 Control:12
## 1st Qu.:-0.38947 1st Qu.:-0.35612 1st Qu.:2021-10-11 Cooling:12
## Median :-0.17933 Median :-0.14491 Median :2021-10-12 Heating:13
## Mean :-0.12610 Mean :-0.08877 Mean :2021-10-13
## 3rd Qu.: 0.09151 3rd Qu.: 0.12617 3rd Qu.:2021-10-18
## Max. : 0.67913 Max. : 0.71289 Max. :2021-10-19
##
## mass_g SVL_mm chamber_time_elapsed_hr date_factor
## Min. : 8.70 Min. :60.00 Min. :0.9333 2021-10-05:6
## 1st Qu.:10.00 1st Qu.:64.00 1st Qu.:1.0000 2021-10-11:8
## Median :11.30 Median :66.00 Median :1.0333 2021-10-12:7
## Mean :11.42 Mean :66.24 Mean :1.1005 2021-10-18:8
## 3rd Qu.:12.60 3rd Qu.:69.00 3rd Qu.:1.1500 2021-10-19:8
## Max. :15.30 Max. :75.00 Max. :1.5333
## length(unique(CEWL_time_series$lizard_ID))## [1] 37NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.
change_LMM_1 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
mass_g + SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_1)## Single term deletions
##
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> 21.929
## treatment 2 42.116
## mass_g 1 20.523
## SVL_mm 1 21.118
## chamber_time_elapsed_hr 1 20.291anova(change_LMM_1)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56693 1.28346 17.8977
## mass_g 1 0.00639 0.00639 0.0891
## SVL_mm 1 0.08371 0.08371 1.1673
## chamber_time_elapsed_hr 1 0.02129 0.02129 0.2969Drop mass first based on low SS and resulting AIC:
change_LMM_2 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_2)## Single term deletions
##
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> 20.523
## treatment 2 40.292
## SVL_mm 1 19.118
## chamber_time_elapsed_hr 1 19.059anova(change_LMM_2)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56389 1.28195 17.5558
## SVL_mm 1 0.02860 0.02860 0.3917
## chamber_time_elapsed_hr 1 0.03428 0.03428 0.4695drop SVL:
change_LMM_3 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_3)## Single term deletions
##
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
## npar AIC
## <none> 19.118
## treatment 2 38.317
## chamber_time_elapsed_hr 1 17.510anova(change_LMM_3)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56421 1.28211 17.8795
## chamber_time_elapsed_hr 1 0.02549 0.02549 0.3554drop chamber_time_elapsed_hr:
change_LMM_4 <- lme4::lmer(data = lizard_coeffs,
rate_change ~
treatment +
(1|date_factor))
drop1(change_LMM_4)## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## rate_change ~ treatment + (1 | date_factor)
## npar AIC
## <none> 17.510
## treatment 2 38.876anova(change_LMM_4)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.5645 1.2823 18.286The simplest model is our best/final model.
Check assumptions:
plot(change_LMM_4)
Assumptions look good.
Save with p-values:
change_LMM_best <- lmerTest::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
(1|date_factor))
summary(change_LMM_best)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rate_change ~ treatment + (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 16.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94436 -0.62201 -0.04392 0.39606 1.77788
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01179 0.1086
## Residual 0.07012 0.2648
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.2093 0.0908 12.1290 -2.305 0.039599 *
## treatmentCooling -0.1793 0.1085 29.9201 -1.653 0.108748
## treatmentHeating 0.4475 0.1068 30.3121 4.191 0.000222 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnC
## tretmntClng -0.593
## tretmntHtng -0.611 0.504anova(change_LMM_best, type = "3", ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 2.5102 1.2551 2 30.708 17.898 7.034e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Treatment group explained a significant amount of the variation (F(3,30)=18.29, p < 0.0001).
Export:
broom.mixed::tidy(change_LMM_best) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/change_CEWL_best_model.csv")also get pairwise differences for change ~ tmt, which I will use to make the boxplot later:
emmeans <- data.frame(emmeans(change_LMM_best, "treatment"))
emmeans## treatment emmean SE df lower.CL upper.CL
## 1 Control -0.2093260 0.09102635 13.09716 -0.40582833 -0.01282375
## 2 Cooling -0.3886322 0.09189944 12.94683 -0.58725181 -0.19001265
## 3 Heating 0.2381211 0.08845166 11.97139 0.04535039 0.43089180emmean_ps <- data.frame(pairs(emmeans(change_LMM_best, "treatment"))) # p values
emmean_ps## contrast estimate SE df t.ratio p.value
## 1 Control - Cooling 0.1793062 0.1087886 30.33842 1.648208 2.414730e-01
## 2 Control - Heating -0.4474471 0.1074827 30.69471 -4.162968 6.708347e-04
## 3 Cooling - Heating -0.6267533 0.1082740 31.09322 -5.788584 6.511073e-06Figures
COLOR
cool <- c(brewer.pal(11, "RdBu")[c(11)])
heat <- c(brewer.pal(11, "RdBu")[c(2)])
control <- c(brewer.pal(9, "Greys")[c(5)])
cntrls <- c(brewer.pal(9, "Greys")[c(7, 5, 3, 5, 4, 5, 3, 6, 7, 3, 6)])cloacal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = calibrated_cloacal_temp,
color = treatment) +
geom_point(size = 0.01,
alpha = 0.2) +
geom_smooth(method = "loess",
se = F,
size = 2,
na.rm = TRUE) +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = "Internal Body Temperature (°C)",
color = "Treatment", se = FALSE) -> ctemp_time## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.ctemp_time## `geom_smooth()` using formula = 'y ~ x'
dorsal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = (calibrated_dorsal_temp),
color = treatment) +
geom_point(size = 0.01,
alpha = 0.2) +
geom_smooth(method = "loess",
se = F,
size = 2,
na.rm = TRUE) +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x= "Time (seconds)" ,
y= "Surface Body Temperature (°C)",
color ="Treatment") -> dtemp_time
dtemp_time## `geom_smooth()` using formula = 'y ~ x'
CEWL (raw) ~ time
ggplot(CEWL_time_series) +
aes(x= time_sec,
y = CEWL_g_m2h,
color = treatment) +
geom_point(size = 0.01,
alpha = 0.25) +
geom_smooth(method = "loess",
size = 2.5,
se = F) +
theme_classic() +
theme(text = element_text(size = 15)) +
theme(axis.text.x = element_text(size = 10)) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = bquote('CEWL (g/'*m^2*'h)'),
color = "Treatment",
se = FALSE) +
scale_x_continuous(limits = c(0, 925),
breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time## `geom_smooth()` using formula = 'y ~ x'
legend
This code uses one of our plots to make a large, detailed legend, which we can save and use in ggarrange later.
ggplot(temp_time_series) +
geom_smooth(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
color = treatment),
method = "lm",
formula = y ~ poly(x, 2),
size = 1,
se = F) +
geom_point(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
fill = treatment,
shape = treatment),
size = 4,
alpha = 1) +
scale_x_continuous(limits = c(14, 39),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = c(0.5,0.5)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") -> fancy_legend
fancy_legend # view
# save
LEGEND <- as_ggplot(get_legend(fancy_legend))
# export
ggsave(filename = "legend_only.pdf",
plot = LEGEND,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 40, height = 30)CEWL (raw) ~ cloacal temperature
ggplot(temp_time_series) +
geom_point(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
shape = treatment,
color = treatment),
size = 0.01,
alpha = 0.1) +
geom_smooth(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
color = treatment),
method = "lm",
formula = y ~ poly(x, 2),
size = 1.5,
se = F) +
scale_x_continuous(limits = c(12.5, 40),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
# annotations for cooling group
annotate(geom = "point", x = 14.3, y = 5,
size = 4, pch = 25, fill = cool) +
annotate("text", x = 12.6, y = 5, label = "T[15]", parse = T,
size = 5, family = "sans", color = cool) +
annotate(geom = "point", x = 35.2, y = 11.8,
size = 4, pch = 25, fill = cool) +
annotate("text", x = 36.7, y = 11.8, label = "T[1]", parse = T,
size = 5, family = "sans", color = cool) +
# annotations for heating group
annotate(geom = "point", x = 15.9, y = 14,
size = 4, pch = 24, fill = heat) +
annotate("text", x = 14.3, y = 14, label = "T[1]", parse = T,
size = 5, family = "sans", color = heat) +
annotate(geom = "point", x = 38.3, y = 25.2,
size = 4, pch = 24, fill = heat) +
annotate("text", x = 40, y = 25.2, label = "T[15]", parse = T,
size = 5, family = "sans", color = heat) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE) -> CEWL_ctemp
CEWL_ctemp
# save
ggsave(filename = "CEWL_int_body_temp.pdf",
plot = CEWL_ctemp,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 150, height = 100)CEWL Control
# get control data only
temp_time_series %>%
dplyr::filter(treatment == "Control") %>%
# pipe into ggplot
ggplot() +
aes(x = calibrated_cloacal_temp,
y = CEWL_g_m2h,
color = lizard_ID) +
geom_point(size = 0.1,
alpha = 0.1) +
geom_smooth(method = "lm",
size = 1,
se = F) +
ylim(0,32) +
xlim(25,31) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = cntrls) +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
color = "Treatment",
se = FALSE) -> CEWL_control
CEWL_control## `geom_smooth()` using formula = 'y ~ x'
# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
plot = CEWL_control,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)## `geom_smooth()` using formula = 'y ~ x'CEWL (relative) ~ time
# check n lizard sample size
CEWL_time_series %>%
dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
group_by(lizard_ID) %>%
summarise(max(relative_CEWL)) %>%
nrow()## [1] 37# plot
ggplot() +
geom_point(data = CEWL_time_series,
aes(x= time_min,
y = relative_CEWL,
shape = treatment,
color = treatment),
size = 0.01,
alpha = 0.1) +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_smooth(data = CEWL_time_series,
aes(x = time_min,
y = relative_CEWL,
color = treatment),
method = "lm",
formula = y ~ poly((x), 2),
size = 1.5,
se = F) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
annotate(geom = "point", x = 15, y = -3.9,
size = 4, pch = 21, fill = control) +
annotate(geom = "point", x = 15, y = -7.5,
size = 4, pch = 25, fill = cool) +
annotate(geom = "point", x = 15, y = 3,
size = 4, pch = 24, fill = heat) +
scale_x_continuous(limits = c(1, 15),
breaks = seq(1, 15, 2)) +
labs(x = "Time (minutes)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
CEWLr_time## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).
# save
ggsave(filename = "CEWL_relative_time.pdf",
plot = CEWLr_time,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 150, height = 100)## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Removed 4 rows containing missing values (`geom_point()`).CEWL Rate of Change
ggplot() +
geom_jitter(data = lizard_coeffs,
aes(x = treatment,
y = rate_change,
shape = treatment,
color = treatment,
fill = treatment),
size = 0.8,
alpha = 0.3,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = emmeans,
aes(x = treatment,
y = emmean,
color = treatment,
group = treatment,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 0.9) +
geom_point(data = emmeans,
aes(x = treatment,
y = emmean,
shape = treatment,
fill = treatment),
color = "black",
size = 4) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
annotate("text", x = 1, y = 0.42, label = "A", size = 3) +
annotate("text", x = 2, y = 0.3, label = "A", size = 3) +
annotate("text", x = 3, y = 0.9, label = "B", size = 3) +
labs(x = "" ,
y = expression(Delta ~ 'CEWL '*min^-1*''),
color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig
# save
ggsave(filename = "CEWL_change_min.pdf",
plot = rate_change_CEWL_fig,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)Experimental Design Skematic
We created reference values to use to make a figure that would show each step of our experiment.
exp_design <- read.csv("./Data/exp design.csv",
header = TRUE,
fileEncoding="UTF-8-BOM")ggplot(exp_design) +
aes(x = time_min,
y = temp_C,
color = treatment) +
geom_line(size = 1) +
geom_vline(xintercept = 60,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(30, 37, family = "sans"), label = "a", color = "black", size = 3) +
geom_vline(xintercept = 65,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(62.5, 37, family = "sans"), label = "b", color = "black", size = 3) +
geom_vline(xintercept = 80,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(72.5, 37, family = "sans"), label = "c", color = "black", size = 3) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.text.x = element_blank(),
#axis.line.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time" ,
y = "Body Temperature (°C)",
color = "Treatment") +
annotate(geom = "point", x = 80, y = 25,
size = 3, pch = 21, fill = control) +
annotate(geom = "point", x = 80, y = 20,
size = 3, pch = 25, fill = cool) +
annotate(geom = "point", x = 80, y = 30,
size = 3, pch = 24, fill = heat) +
scale_y_continuous(limits = c(14, 37),
breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig
# save
ggsave(filename = "methods_schematic.pdf",
plot = methods_schmatic_fig,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)Save Grouped Figures
Relative CEWL & Rate of Change
multi_fig_relative <- ggarrange(CEWLr_time,
ggarrange(rate_change_CEWL_fig, LEGEND,
nrow = 1, ncol = 2,
widths = c(2,1),
align = "hv"
),
nrow = 2, ncol = 1,
heights = c(2,1),
labels = c("A", "B"))## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).multi_fig_relative
ggsave(filename = "CEWL_time_multi_fig.pdf",
plot = multi_fig_relative,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 120, height = 160)Raw CEWL ~ Cloacal temp & control zoom-in
multi_fig_CEWL_by_temp <- ggarrange(CEWL_ctemp,
ggarrange(CEWL_control, LEGEND,
nrow = 1, ncol = 2,
widths = c(2, 1),
align = "hv"
),
nrow = 2, ncol = 1,
heights = c(2,1),
labels = c("A", "B"))## `geom_smooth()` using formula = 'y ~ x'multi_fig_CEWL_by_temp
ggsave(filename = "CEWL_temp_multi_fig.pdf",
plot = multi_fig_CEWL_by_temp,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 120, height = 160)